Minimal Mechanisms of Microtubule Length Regulation in Living Cells

The microtubule cytoskeleton is responsible for sustained, long-range intracellular transport of mRNAs, proteins, and organelles in neurons. Neuronal microtubules must be stable enough to ensure reliable transport, but they also undergo dynamic instability, as their plus and minus ends continuously switch between growth and shrinking. This process allows for continuous rebuilding of the cytoskeleton and for flexibility in injury settings. Motivated by in vivo experimental data on microtubule behavior in Drosophila neurons, we propose a mathematical model of dendritic microtubule dynamics, with a focus on understanding microtubule length, velocity, and state-duration distributions. We find that limitations on microtubule growth phases are needed for realistic dynamics, but the type of limiting mechanism leads to qualitatively different responses to plausible experimental perturbations. We therefore propose and investigate two minimally-complex length-limiting factors: limitation due to resource (tubulin) constraints and limitation due to catastrophe of large-length microtubules. We combine simulations of a detailed stochastic model with steady-state analysis of a mean-field ordinary differential equations model to map out qualitatively distinct parameter regimes. This provides a basis for predicting changes in microtubule dynamics, tubulin allocation, and the turnover rate of tubulin within microtubules in different experimental environments.


Introduction
Microtubules (MTs) are protein filaments that provide structure to cells and allow them to take on diverse shapes.For instance, neurons are long-range cells where protein transport must occur across long distances, and this is crucially supported by the microtubule cytoskeleton (Kapitein andHoogenraad, 2011, 2015;Kelliher et al, 2019;Rolls et al, 2021).These polymers must therefore be stable enough to support reliable transport of various proteins and organelles.However, microtubules also have dynamic ends, which contribute to the continuous rebuilding of the cytoskeleton.For example, turnover experiments in segments of dendrites of Drosophila neurons show that microtubule polymers, comprised of tubulin, turn over every few hours in these cells (Tao et al, 2016).This is due to the switch-like growth behavior that both ends exhibit, where periods of growth transition suddenly into shrinking behavior (catastrophe) and back to growth (rescue).This behavior, called dynamic instability (Mitchison and Kirschner, 1984), is especially important when neurons are confronted with stress and injury (Rolls et al, 2021).
Given the complex and intrinsically stochastic dynamics of the microtubule cytoskeleton, many prior mathematical models of MT dynamics have been developed to describe it.These works have used various mathematical frameworks, but usually their primary goal was to reproduce fine details of dynamic instability of MTs observed from in vitro experiments.The approaches proposed included two-phase models of microtubule dynamics, where polymers can switch between growing and shrinking phases and the dynamics are primarily described using coupled partial differential equations (Hill, 1984;Dogterom and Leibler, 1993;Govindan and Spillman, 2004;Akhmanova and Kapitein, 2022;Honoré et al, 2019;Barlukova et al, 2018).Follow-up studies proposed stochastic and continuous models of microtubule growth and catastrophe incorporating GTP hydrolysis mechanisms (Flyvbjerg et al, 1996;Deymier et al, 2005;Margolin et al, 2006;Antal et al, 2007;Hinow et al, 2009;Mazilu et al, 2010;Buxton et al, 2010;Padinhateeri et al, 2012;Barlukova et al, 2018) and reflecting the observation that catastrophe may be a multi-step process (Bowne-Anderson et al, 2013).More complex computational models combined known or assumed mechanical, thermodynamic, and kinetic mechanisms for tubulin assembly and bending, GTP hydrolysis, or molecular dynamics of tubulin-tubulin interactions to understand the structure and fluctuation of microtubule tips during dynamic instability (Molodtsov et al, 2005;Ji and Feng, 2011;Margolin et al, 2012;Zakharov et al, 2015).
Despite the dynamic instability of individual filaments, microtubule networks are remarkably stable and long-lived, particularly in neurons (Rolls et al, 2021).One population-scale quantity that has attracted attention from mathematical modelers is the distribution of microtubule lengths.While some studies observed and characterized the transition between parameter regimes yielding unbounded growth and finite steady-state MT lengths, fewer models incorporated limited tubulin availability or the impact of tubulin diffusion on MT length control (Dogterom et al, 1995;Deymier et al, 2005;Hinow et al, 2009;Margolin et al, 2006;Buxton et al, 2010).Even fewer considered mechanisms where MT regulation may arise from the dependence of certain parameters on the length of the MT.For instance, the model proposed in Mazilu et al (2010) assumed that the attachment and detachment rates of GTP monomer units to the filament depend linearly on the MT length, motivated by the physical constraint that MT lengths cannot be negative, nor exceed a maximum length.In Rank et al (2018), the authors included a tubulin resource limitation in their model by allowing the filament polymerization rate to decrease linearly with increasing MT length, while in Bolterauer et al (1999), they assumed that the catastrophe probability is an increasing function of length.However, a comprehensive characterization of how multiple mechanisms such as competition for tubulin resources and length-dependent catastrophe of MTs contribute to MT regulation is lacking.
In addition, a few of these previous works considered realistic living cell behavior, such as regulation of MT length in the mitotic spindle (Dogterom et al, 1995;Rank et al, 2018) or the role of cell geometry, nucleation center, and cell edges on dynamic instability (Deymier et al, 2005;Margolin et al, 2006;Buxton et al, 2010).These studies focused on models of MTs growing from a centrosome or nucleation center and thus exhibiting dynamics at the plus end only.Most studies used kinetic parameters for MT dynamics extracted from in vitro studies or otherwise simulated the MT behavior for large parameter ranges.However, arrays of non-centrosomal MTs are key for the specialized functions of many cells (Keating and Borisy, 1999; Bartolini and Gundersen, 2006;Sanchez and Feldman, 2017;Akhmanova and Kapitein, 2022).In many non-centrosomal arrays, minus ends are localized to specific regions by nucleation or capping proteins including the g-tubulin ring complex and CAMSAPs, while the plus end remains free and dynamic.In neurons, MTs form an overlapping array with plus and minus ends staggered throughout axons and dendrites (Baas and Lin, 2011;Kapitein and Hoogenraad, 2015).Plus ends grow and shrink throughout axons and dendrites (Baas et al, 2016;Rolls et al, 2021), and at least in Drosophila fruit flies and in zebrafish, minus ends also exhibit dynamic instability (Feng et al, 2019).
Addressing the question of how MT length is regulated in a living cell is especially important when studying cytoskeleton dynamics during cell growth or response to injury.While previous complex models of microtubule behavior have provided useful insights for the understanding of dynamic instability in vitro, studying wholecell cytoskeleton dynamics may not require the same level of small-scale biochemical details.Understanding how local fluctuations in MTs impact the overall stability of the MT cytoskeleton in a healthy living cell requires a minimal modeling approach with parameters that can be validated with experiments.An advantage of a minimal modeling strategy is that it would not attempt to incorporate the effects of individual molecules on microtubule growth at tips, which is still actively studied experimentally.Such a modeling framework would also enable further investigation of microtubule dynamics during developmental or injury processes.
Even with such a minimalistic modeling framework, there are numerous challenges.For instance, it is difficult to fully observe data on global microtubule behavior (such as the polarity, density, or length of microtubules) in living cells.Local information on microtubule behavior can be more readily quantified, but may only reflect partial observations.For example, run lengths of microtubule ends may be cut off due to experimental time windows, and parameters extracted from such measurements are likely to represent effective quantities, which already incorporate mechanisms of MT length regulation in an unknown way.
In summary, the main goal of this work is to investigate two minimal lengthlimiting factors-limitation due to resource (tubulin) constraints and limitation due to catastrophe of large-length microtubules-and to understand how these different constraints then shape steady-state statistics of microtubule dynamics.Moreover, we study how steady-state length distributions evolve in response to parameter changes.Since the action of these limiting factors is difficult to capture directly, a predictive stochastic model enables the study of a variety of experimentally observable emergent properties.The stochastic modeling framework has the advantage that it can capture fluctuations in behavior within and among different microtubule trajectories, however it requires many parameters in a parameter space that is difficult and computationally costly to navigate.We therefore propose a coarse-grained deterministic (ordinary differential equations) model that explicitly addresses certain stochastic effects in a novel way.We find that the average MT length is not sufficient to identify all parameters of interest in the stochastic model.However, the deterministic framework helps restrict our study to parameter combinations that generate target MT lengths, thus reducing the feasible parameter space.As a consequence, we can systematically investigate the large parameter space and make faithful predictions in multiple parameter regimes.Informing the stochastic model in this way allows us to distinguish between the contributions of different mechanisms of microtubule length regulation, and to provide predictions for experiments that may soon be possible, such as controlling tubulin concentration in vivo.

Overview of modeling decisions
We propose a stochastic model that tracks the dynamics of several microtubules along one spatial dimension.This one-dimensional space assumption is appropriate for microtubules extending along a segment of an axon or dendrite of a neuron.We model polymerization and depolymerization dynamics at both the plus and minus end of the microtubules.Incorporating the dynamics at both ends is seldom studied in prior models, especially those motivated by mitotic spindle dynamics (Deymier et al, 2005;Dogterom et al, 1995;Rank et al, 2018;Bowne-Anderson et al, 2013).We are motivated by systems in which microtubules are dynamic at both ends, consistent with recent evidence on dynamic instability at both microtubule ends in neurons of Drosophila fruit flies and zebrafish (Feng et al, 2019;Thyagarajan et al, 2022).In addition to overall growth and shrinking of microtubules, these systems may also exhibit treadmilling (growth at one end and shortening at the other end) (Rodionov and Borisy, 1997;Margolis and Wilson, 1998;White et al, 2015).As motivated above, we are interested in understanding how MT length distributions are regulated by one or two qualitatively different growth control mechanisms: a regime that is limited by competition for tubulin, and a regime that is limited by intrinsic length-driven limitation of catastrophe.
The assumption that microtubules compete for tubulin resources is widely accepted and has been implemented in previous mathematical models (Dogterom et al, 1995;Deymier et al, 2005;Hinow et al, 2009;Margolin et al, 2006;Buxton et al, 2010); in our framework, this potential limitation due to tubulin availability impacts the growth velocity of the microtubules at both ends.While we do not include explicit spatial diffusion of tubulin, we consider separate pools of free and unavailable tubulin.Unlike previous models that included a linear dependence of polymerization rates on free tubulin concentration (Bolterauer et al, 1999;Rank et al, 2018), we use Michaelis-Menten kinetics to model a smooth dependence of the growth velocity on the free tubulin pool.
Modeling a length-dependent catastrophe mechanism is motivated by observations from in vitro and in vivo data (Gardner et al, 2013).Experiments have shown that catastrophe is suppressed in short newly polymerizing microtubules (Gardner et al, 2011), and the initial simple hypotheses of stochastic changes in GTP cap length leading to catastrophe could not account for the increase in catastrophe frequency as microtubules grow.Two types of mechanisms have been proposed to account for the increased catastrophe frequency as microtubules grow.The first is the "antenna" mechanism for kinesin-stimulated catastrophe, which hypothesizes that a depolymerizing kinesin like kinesin-8 drops down on a microtubule and walks to the plus end where it removes subunits to promote catastrophe.Longer microtubules would receive more motors and thus have more subunit removal at the plus end (Varga et al, 2006(Varga et al, , 2009)).The second major hypothesized mechanism posits that, as microtubules grow, the plus ends become more ragged, with some protofilaments lagging behind while others keep extending.This accumulation of more variable protofilament lengths as microtubules elongate is supported by observations on plus ends in vitro and in vivo (Coombes et al, 2013), and different modeling approaches support the idea that this could be responsible for length-dependent catastrophe (Alexandrova et al, 2022;Coombes et al, 2013).Here, we model length-dependent catastrophe of microtubules by imposing a rate of switching from growth to shrinking that depends linearly on microtubule length.
Mathematical models have been developed to probe these and other mechanisms of catastrophe and rescue in microtubule dynamics.Many of these either involved a high level of molecular detail (Flyvbjerg et al, 1996;Deymier et al, 2005;Margolin et al, 2006;Antal et al, 2007;Hinow et al, 2009;Mazilu et al, 2010;Buxton et al, 2010;Padinhateeri et al, 2012;Barlukova et al, 2018;Bowne-Anderson et al, 2013;Molodtsov et al, 2005;Ji and Feng, 2011;Margolin et al, 2012;Coombes et al, 2013;Zakharov et al, 2015), or incorporated complexity through memory effects, such as age-dependent microtubule catastrophe (Barlukova et al, 2018) or location-dependent rescue (Fees and Moore, 2019).However, including this level of detail or memory comes with the cost of analytical intractability.Since our goal is to investigate largescale microtubule behavior in whole cells, and in vivo experimental data at this scale is limited, we proceed with memory-less versions of the switching mechanisms, and adopt a continuous-time Markov chain (CTMC) framework, which is summarized by Figure 1.This enables a detailed understanding of qualitative behaviors arising from different parameter regimes and an ability to predict the impact of realistic experimental perturbations.

Sh rinking Growth
Free tubulin

Unavailable tubulin
Microtubule tubulin pool Microtubule end states < l a t e x i t s h a 1 _ b a s e 6 4 = " Y Q q 7 o a q P 9 u O p y H j B 2 Z 4 5 l a t e x i t s h a 1 _ b a s e 6 4 = " 2 P Q 4 0 q p y m t w m H 3 e o 9 x 1 a V h 9 l h + 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " s 9 G 2 u f D q h P f p z r 4 N G D D y E L P M Z 0 E = " > A A A B + H i c b V D L S g N B E J z 1 G e M j q x 6 9 D A b B i 2 F X g n o M e v E Y w T w g W Z f Z y W w y Z P b B T I 8 Q l / 0 S L x 4 U 8 e q n e P N v n C R 7 0 M S C h q K q m + 6 u I B V c g e N 8 W y u r a + s b m 6 W t 8 v b O 7 l 7 F 3 j 9 o q 0 R L y < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 C h s 2 y Z F p q s t x l w F 4 M / x J 1 6 q g G c = " In growth, the MTs utilize the available tubulin pool; while shrinking, they replenish tubulin into the unavailable pool.

Overview of the paper
In Figure 2, we show a schematic that outlines our parameterization processes, involving both a stochastic CTMC model and a coarse-grained ordinary differential equations (ODE) model.Section 2 outlines the CTMC model that describes polymerization and depolymerization of microtubules under various length regulation constraints in a segment of a dendrite.The stochastic model easily relates to experiments, which include observables like minus-end trajectories, average MT length, and growth speed distributions.However, it is difficult to understand the qualitative behavior of the MTs given the large parameter space of this model.Therefore, Section 3 introduces the reduced ODE model of MT growth dynamics that we analyze in order to propose a strategy for parameterization of the stochastic model and to identify relevant different parameter regimes.We treat some experimental parameters as more reliably-observed than others, for example, the growth and shrinkage speeds at both MT ends.We use experimental data in Drosophila dendrites as our primary guide for parameterization.Then, in order to target certain desired emergent properties like average MT length, we derive a protocol for assigning values to unobservable parameters in Section 3.3.Given this protocol, we then investigate observables for prediction like fluorescence turnover (see Figure 2) in Section 4 based on the parameterized stochastic model in qualitatively different regimes of the MT behavior and under potential experimental manipulations.We discuss our results in Section 5.

Stochastic model of microtubule dynamics
To investigate the dynamics of growing and shrinking microtubules, we use a stochastic CTMC model.In each simulation, we assume a fixed number N of microtubules.This is similar to previous investigations (Margolin et al, 2006;Jonasson et al, 2020), and consistent with the idea that microtubule populations are tightly regulated in vivo.

Observables for predic�on
Sec�ons 2.1, 4  We consider a one-dimensional domain, modeling a neuronal dendrite segment that has uniform polarity.This means that each microtubule is summarized by the pair (x + (t), x − (t)), where x + (t) and x − (t) refer to the position of the plus and minus end, respectively.Note that we assume x + (t) ≥ x − (t) for all t.Initially, each MT end is randomly assigned to be in one of two states: growth, where the microtubule end is polymerizing, and shrinking, where the microtubule is depolymerizing.Each microtubule end can switch from one state to another at some rate λ.We denote the switching rate from growth to shrinking at the plus end of the MT as λ + g→s , and the rate of switching from shrinking to growth as λ + s→g .In the shrinking state, the MT plus ends shrink at velocity v + s , and in the growth state, the plus ends grow with velocity or v + g .We use similar notation for describing the dynamics at the minus end.A cartoon depiction of our CTMC stochastic model framework is shown in Figure 1.
We consider a small time step ∆t = 1 second for the stochastic model implementation.We then make the assumption that there can be at most one shrinking event in one time step.All possible events that can occur in a time step are shown in Figure A2 in Appendix A. Depending on the state of the MT end and the time spent in each state in that time step ∆t, we implement shrinking events and then growth events at the MT ends.The MT end positions are then updated according to these events, and we proceed to the next time step.Additional details on the implementation of the stochastic model are provided in Appendix A.
In vivo experimental measurements for plus and minus end growth velocities and run times have been obtained for microtubules in Drosophila dendrites in Feng et al (2019), while run times and velocities in shrinking states have also been investigated in this system (data unpublished).We can use these experimental quantities to parameterize our CTMC model, allowing us to simulate MT dynamics and observe how filament behavior emerges over time.Figure 3a shows trajectory paths for the plus and minus ends of one microtubule (left) and lengths for 20 microtubules (right) over the course of five hours of simulation time based on the available experimental measurements in Drosophila dendrites.The mean microtubule length over time is shown in red, and the red cloud illustrates the interquartile range of lengths over time.The left panel shows that the plus end is moving further to the right through time, while in the right panel, each MT length (grey) is increasing overall throughout time.The length of the dendrites in some Drosophila neurons are on the order of tens of microns (Stone et al, 2010), so microtubule lengths on the order of 1000 µm are not physiologically reasonable for this system.Imposing a domain restriction in this model would lead to microtubules extending throughout the length of the dendrites, which is also unreasonable for this application.
We thus propose two mechanisms to regulate unbounded MT growth: tubulindependent growth and length-dependent catastrophe.Figure 1 illustrates how the two mechanisms are incorporated into our CTMC simulation, where growth velocity is a function of free tubulin and the switching from growth to shrinking is dependent on the length of a microtubule (L).Specifically, we assume that the switching rate from growth to shrinking satisfies the following conditions.Definition 2.1.A length-dependent catastrophe rate is defined in terms of four quantities: L 0 > 0, the characteristic MT length; λ ± g→s > 0, the plus-end (respectively, minus-end) growth-to-shrinking phase switch rate when a MT has length L 0 ; λ min ∈ (0, λ ± g→s ), the smallest possible value of the catastrophe rate; γ ≥ 0, a parameter that captures the magnitude of length-dependence in the switch rate; and ϕ : R + → R + , a dimensionless shape function for the length-dependence.
We say that a function ϕ is an admissible length-dependent catastrophe shape function if it is continuous, non-decreasing, and satisfies ϕ(0) = −1, ϕ(1) = 0, and ϕ ′ (1) = 1. (1) When a MT has length L, In this way, the catastrophe rate is always strictly positive and equals λ ± g→s when L = L 0 .Moreover, the derivative of the rate equals γ at the characteristic length L 0 .In practice, we will use ϕ(x) = x − 1 for the shape of the rate function.Note that since ϕ(0) = −1, when γ < (λ ± g→s − λ min )/L 0 we have that the rate exceeds λ min for all L. Otherwise, the minimum rate λ min will be imposed for short MTs.
When including these mechanisms, microtubule growth is facilitated by the availability of tubulin dimers.Tubulin protein can thus either be found in microtubules or be available throughout the neuron for growth.We assume that growth is dependent on the amount of free tubulin (F ), and that a shrinking event releases tubulin and thus leads to an increase in F .To account for possible slow diffusion of tubulin throughout the length of the dendrite, we assume that after shrinking, the newly-released tubulin is unavailable for a short period of time, which we denote as unavailable tubulin, U .After some time τ tub , the unavailable tubulin U is made available and becomes free tubulin, F .As stated in section 1.1, growth is modeled using Michaelis-Menten kinetics to yield velocities that saturate in rich tubulin environments and that decrease smoothly to a small growth velocity when free tubulin F is low.

Stochastic model observables
Our goal is to use the stochastic model framework outlined in this section to simulate MT dynamics that achieves target measurements such as mean growth speed, mean growth segment duration, and mean MT length (Figure 2).Besides these quantities, the model proposed also predicts several emergent measurements that we discuss here.Some of these measurements are quantities that we predict and are not currently available from experiments, while others can be validated with existing experimental work (Figure 2).
For instance, the model quantifies the allocation of tubulin resources into available, unavailable, and microtubule pools.These different tubulin pools are currently very challenging to observe experimentally.In addition, the model can capture a predicted movement direction and speed of the MTs, which emerge as a result of the overall polymerization and depolymerization dynamics at both ends.Figure 3b shows the predicted trajectory of a single MT in the minus end direction, which suggests some treadmilling of the microtubule.
The framework proposed can also be used to model and validate the dynamics using photoconversion experiments, which represent an approach for measuring microtubule stability and turnover in vivo.
In axons and dendrites of neurons, photoconversion experiments tag tubulin in a region of interest and observe the fluorescence in that region at subsequent time points (Rolls et al, 2021) (see Figure 4).The fluorescence measured in the window corresponds to tagged tubulin protein that is trapped in a portion of the MT that remains polymerized.Faster decay of the fluorescence in the region of interest indicates more dynamic MTs, which depolymerize on short timescales (Rolls et al, 2021).
We can extract this measurement of MT turnover using our stochastic model as follows.We first select a photoconversion time (t = 120 minutes), which denotes the time in the simulation when tubulin is tagged in a 10 µm window of interest.We then calculate the microtubule content distribution throughout space at the photoconversion time and choose placement of the window of interest in a way that maximizes the tubulin content (from MTs) in the window.This choice allows us to observe the dynamics of as many microtubule segments as possible in the photoconversion window.At all subsequent times, we track the location of the MT segments that were tagged during the simulated photoconversion and determine the portions of these segments that remain within the window and have not yet depolymerized.As in experiments, we are interested in the relative fluorescence intensity remaining in the window of interest as a function of time (Tao et al, 2016;Rolls et al, 2021), and therefore normalize the fluorescence measurement at the photoconversion time to 1 in all simulation trials.
For the rest of the paper, we seek to understand what mechanisms affect microtubule behavior, which can be quantified through our stochastic model observables.However, the stochastic simulation relies on many parameters related to growth and shrinking dynamics, and how those parameters depend on each other and influence microtubule growth dynamics is unknown.Additionally, running CTMC simulations with a small timestep is computationally expensive.In the next section, we explore a continuous model of tubulin allocation and microtubule growth to understand how model observables change with different stochastic parameters.Since we ultimately wish to qualitatively match experimental data such as microtubule growth speed and growth duration, we would like to be able to tune our stochastic parameters to match these data.
3 Reduced deterministic model of microtubule dynamics

Statement of ODE model and justification
In order to identify and characterize fundamentally distinct parameter regimes, and to establish a protocol for fixing parameters for our investigation, we introduce a deterministic ODE model with a reduced version of the dynamics.Our deterministic approximation is a compartmental model, using the flux of tubulin in and out of the microtubule pool to define the terms of the ODE.To simplify the analysis, we will neglect the unavailable pool of tubulin.In this way, we need only track the amount of tubulin in microtubules at time t, M (t), with the amount of free tubulin satisfying , where T tot is the (constant) total amount of tubulin in the system.Meanwhile, let G + (t) and G − (t) denote the number of MTs that are in growth phase at the plus end and minus end, respectively.We assume a fixed number of MTs, N , and let N denote the average MT length.Similarly, we let g N correspond to the fraction of MTs in growth phase at each end.Many of the features of the deterministic model are a direct translation or slight generalization of what we use for stochastic model components.One example is the form of the length-dependent catastrophe, which is here expressed in terms of a dimensionless shape function ϕ(x) (see Definition 1).On the other hand, there are aspects of the mean-field model that are not directly represented in the stochastic version.For example, we introduce a function ψ α (x) that phenomenologically captures the tendency of microtubules to hit zero-length during shrinkage runs, using the shape parameter α to capture the effect for different MT length distributions.We express our ODE model using qualitative properties of ϕ and ψ α , but a specific definition for ψ α can be found in Equation 11 at the end of our model derivation.Definition 3.1.Given a shape parameter α ≥ 1, let the active MT proportion function ψ α : R + → R + be an increasing and continuously differentiable function satisfying (3) Definition 3.2.Let ϕ and ψ α satisfy the conditions of Definitions 2.1 and 3.1, respectively.We define our mean-field approximation of the stochastic model to be the system of ODEs: (4) The model parameters are described in Table 1.By dividing the equations in model ( 4) by the constant number of MTs N , we obtain the following equations that we focus our subsequent analysis on: (5) In Equations ( 5), the fraction of growing microtubules on either the plus or minus end (denoted as g + (t) and g − (t), respectively) is dictated by the flux of microtubules switching into and out of the growing state.The switch from shrinking to growth on either end is given by rate λ ± s→g , which we assume to be constant.The switch from growth to shrinking is driven by the length-dependent catastrophe mechanism in Definition 2.1, where the rate of catastrophe is linearly dependent on the MT length L(t), with a minimum rate λ min imposed for short MTs.
The flux of tubulin into the microtubule pool in Equations ( 5) is governed by the growth rate of microtubules at the plus and minus ends.We denote the maximum growth speeds by v + g and v − g .However, the actual growth rate can be tempered by the scarcity of available tubulin, which we model with the Michaelis-Menten term in Equations ( 5), governed by the half-rate parameter F 1/2 .In Equations ( 4), before the normalization by the number of MTs N , this corresponds to: Tubulin flux into MT pool: The flux out of the MT pool is difficult to directly model.We assume that in reality the shrinkage rate has a constant mean v − g and that this is maintained until the MT hits zero length.In our model, the "shrinking state" continues until a switch to growth state occurs (all at Markovian rate λ ± s→g ).The deterministic model is an average over all microtubules, which is to say a weighted average of MTs with nonzero length (with shrinking velocity v ± s ) and of MTs with zero length (hence with shrinking velocity zero).We thus need a term that models the reduction in overall tubulin loss when some MTs in shrinking state have length zero.
Indeed, preliminary simulation work showed that when there is no MT length dependence of the catastrophe rate, the MT length distribution looks exponentially distributed (Figure 6), and certainly strictly decreasing with length.On the other hand, when MT-length dependence is strong, the distribution is no longer monotone: with a mode that is closer to the median, but is nevertheless quite wide.Throughout our simulations, when the average MT is within a few multiples of the average run length of a shrinking state, then we expect simulations to contain a non-trivial number of shrinking-state MTs with zero length.
To capture this translation from an individual-based stochastic model to a meanfield ODE dynamical system, we introduce a family of dimensionless functions ψ α : R + → [0, 1], indexed by a shape parameter α > 0, that we use to approximate the proportion of MTs that have non-zero length when the mean of the length distribution is L. Let L crit be a characteristic length, chosen below, at which a microtubule is "at risk" of becoming zero length.Then if Y is a random length drawn from the MT-length distribution with shape α and mean L, we write On the right-hand side, we see that this is the probabilistic complement of the cumulative distribution function (cdf) of Y .When the distribution of Y is a member of a shape-and-scale family, we can express its cdf in terms of the cdf of a unit-scale member of the family that has the same shape parameter.Suppose that X has this unit-scale distribution, and write F α (x) = P α (X ≤ x).Commonly the unit-scale distribution will have a mean that is different than one, so we define µ α = E α (X).Then we can define Y = LX µα and observe that In light of ( 7) this means that we can express ψ α in terms of the unit-scale cdf model: where x = L Lcrit .Note that, by the properties of cumulative distribution functions for strictly positive random variables, the conditions of Definition 3.2 for ψ α will commonly be met.
It remains to choose L crit and a family of probability distributions that emulate the qualitative behavior of the MT-length distributions.Since L crit represents a length for which the MT is "at risk" of being zero length, we choose to let L crit be the length for which there is a 50% chance that a MT will hit zero length before the end of a shrinking run.Since the plus-end runs are exponentially distributed with mean v + s /λ + s→g , an appropriate choice for the critical length is the median of that distribution, so As discussed above, preliminary numerical investigation revealed that when there is no length-dependence in the catastrophe rate, the distribution of MT lengths is qualitatively exponential.As length-dependence becomes more prominent, the distribution is no longer monotone decreasing, but retains very wide support.A natural generalization of the Exponential family is the Weibull family of distributions.This family was originally introduced to model hazard rates (in our context, state-switch rates) that are increasing with time (in our context, with MT length).Moreover, the complement to the cdf has a very simple expression: 1 − F α (x) = exp(−x α ).Since the mean of a Weibull(α, 1) random variable is µ α = Γ(1 + 1 α ), we have the following final definition for ψ α : (11)

Non-dimensionalization and basic mathematical properties
There are several natural candidates for time-and length-scales to use for the purpose of nondimensionalization.A few important length scales are the average lengths of growth and shrinkage runs (the ratio v + g /λ + g→s is the average run length of a plusend growth phase, for example); the total amount of tubulin that is available per microtubule (T tot /N ); and the characteristic length L 0 .We choose the latter, in part because then all quantities can be understood as a fraction of the characteristic MT length scale; then the nondimensionalized dynamics for MT length will be restricted to the interval [0, T tot /(N L 0 )].For the time scale, we choose the average run duration of a plus-end shrinking run (1/λ + s→g ).As we have seen, this quantity plays an important role in establishing the critical length scale of the mean-field approximation for the number of MTs that have non-zero length.
Table 2: Dimensionless groups that appear in the non-dimensional ODE system (12).
Given the choices of length and time scale presented in Table 2, with t = t t char and ℓ = L L char , we have the following nondimensional system: ) where ′ denotes the derivative with respect to the nondimensional time t.Proposition 3.3.Suppose that η ± ≥ 0, T > N ℓ 0 , and all other nondimensional constants in Table 2 are positive.Let O = (0, T /N ) × (0, 1) × (0, 1).Then for any initial condition (ℓ, g + , g − ) in O there exists a unique global solution to the system (12).Moreover, there exists a unique steady-state solution (ℓ * , g + * , g − * ) in the region O. Remark 3.4.We do not have a proof that the steady state is asymptotically stable, or that chaos can be excluded from this three-dimensional system.However, numerical experimentation strongly indicates there is convergence to the unique steady-state.
Proof.First we observe that the unit cube (the domain closure, O) is a forwardinvariant domain for the dynamics.That is to say, the flow from all points on the boundary are to the interior.For example, if either g + = 0 or g − = 0, then their derivatives are 1 and 1 λs→g , respectively, and the flow is into the interior of the domain.On the other hand, suppose that either g + = 1 or g − = 1, the positivity of λ min assures that their derivatives are negative.Meanwhile, 0 < g + , g − < 1 and ℓ = 0 results in ℓ ′ > 0. Similarly, ℓ = T /N results in ℓ ′ < 0. To summarize, if the system takes its value in the interior or edge of any faces of the parallelepiped domain, the flow is to the interior of the domain.Finally, we note that none of the corners are fixed points, since the local dynamics pushes the system to a domain face, from which it enters the interior of the domain.For example, at (0, 0, 0), we have that ℓ ′ = 0, but (g + ) ′ > 0 and (g − ) ′ > 0. The domain is therefore a forward invariant set.Because all functions in the system are continuous in the domain with bounded derivatives, they are Lipschitz continuous.Global existence and uniqueness of solutions follows.
To establish the existence and uniqueness of a steady state, note first that by setting derivatives equal to zero, we can immediately solve for g ± * , if ℓ * is given.In fact, for any value ℓ, the growth versus shrinking phase fraction is in balance at the values g ± * (ℓ) where These values lie strictly between 0 and 1.Moreover, since ϕ(ℓ) is non-decreasing in ℓ, both g + * and g − * are decreasing in ℓ.Turning our attention to the equation for ℓ * itself, the steady-state value must satisfy where Since ψ α (0) = 0, h 0 (0) = 0.As ℓ increases, we note that h 0 (ℓ) is a product of positive increasing functions, and is hence an increasing function itself.Meanwhile, h 1 (0) > 0 and h 1 ( T /N ) = 0 while decreasing in ℓ in between.It follows that h(l) is decreasing with h(0) > 0 and h( T /N ) < 0. By the Intermediate Value Theorem, there exists a unique ℓ * ∈ (0, T /N ) such that h(ℓ * ) = 0, so that ( 14) holds.It follows that (ℓ * , g + * (ℓ * ), g − * (ℓ * )) is a unique fixed point of the system (12).We display a visualization of the fixed point argument in Figure 5 through modified versions of the functions h 0 and h 1 .Since the size of the growth phase population is never zero after t > 0, we can divide Equation 14 through by g + * (ℓ) + g − * (ℓ) vg and define the functions The benefit of this presentation is that the functions H 1 and H 0 split the dependence on the two major parameters we will perturb in our analysis.The function H 1 depends on the total tubulin in the system T tot through T , but not on the magnitude γ of length-dependence in the catastrophe rate.On the other hand, since L char is chosen to be independent of the total tubulin, then H 0 depends on γ (in non-dimensional form, on η) through the g ± * functions, but is not affected by changes in T tot .Therefore, as we perturb these two parameters, they lead to a shift in one function or the other, which we can visualize as in the left panel of Figure 5.For each combination of T tot and γ, there is a unique intersection of H 0 and H 1 , which is known to exist due to the arguments concerning h 0 and h 1 in the proof of Proposition 3.3.16).We use the baseline parameter set in both cases, with the modifications that H 0 is evaluated over a range of γ (blue, red, and green), while the function H 1 is evaluated over a range of total tubulin (light gray to black as T tot increases.)The steady state values (filled circles) occur at intersections of H 0 and H 1 .(Right) The steady state L * is displayed as a function of T tot for the fixed values of γ.The displayed steady state values (filled circles) correspond to those displayed on the left.

Strategy for parameterization
One of our fundamental goals is to be able to observe a system, find a best-fit set of parameters, and then perturb selected model parameters to predict how the system will respond if such a biological experiment was carried out.An example of such an exercise can be seen in the right panel of Figure 5. First, using the procedure described below, we selected parameters in such a way that the "experimentallyinformed quantities" listed in Table 1 were respected while setting L * = 35µm when N = 20 and T tot = 1000µm for each choice of γ ∈ {0, 0.005, 0.03}.The success of the presented procedure can be seen in the left panel of Figure 5, because the functions H 0 and H 1 all intersect at the same location for this fixed value of T tot .Fixing all parameters except for T tot , we then solved for L * as we varied T tot from 750 to 4000 µm.We display the relationship between L * and T tot for three choices of γ.When γ = 0, there is no length-dependence in the catastrophe rate, and therefore no penalty in MTs growing larger.The steady-state average MT length value increases linearly without bound.On the other hand, strong length-dependence in the catastrophe rate results in very little increase in L * when more tubulin is available in the system.Because the values γ = 0, 0.005, 0.03 yield qualitatively different responses in L * to T tot , we adopt these as our standard choices when conducting numerical experiments for the stochastic model in Section 4.
In this section, we outline a strategy for selecting model parameters for the stochastic system relying on quantities that are most directly observable.For parameters that are not known, we can derive quantities that will result in fixed points of the ODE that generate the desired values (see Figure 2).While not a perfect match to the full stochastic model, we have found very good agreement when the total tubulin in the system T tot exceeds the amount required for the MTs to reach their targeted average length.
Referencing the "experimentally-informed quantities" and the "prescribed model parameters" in Table 1, we select our parameters so that the following constraints are met: (C1) the resulting average microtubule length in steady state is L * = L; (C2) the average speed of growth phase polymerization is a fraction of the maximum polymerization speeds in a ratio consistent with the observable quantities v ± g and v ±,max g ; and, (C3) the growth to shrinking phase switch rate at steady-state matches the inverse of the observed growth phase durations τ ± g .Constraints (C2) and (C3) address speeds and durations of MTs in growth phase, however we do not articulate similar constraints for shrinking states.For shrinking state speed, this is because we observed, in preliminary data, that the variance of shrinking phase speeds was significantly less than that of growth phase speeds.We hypothesize that this is because growth phase polymerization is sensitive to tubulin availability, while shrinking phase depolymerization is not subject to an external constraint.We can therefore simply set the model parameters v ± s equal to the experimentally observable parameters v ± s .On the other hand, there is not a directly measureable analogue of the switch rate from shrinking to growth.This is because we experimentally measure the time spent in the shrinking state, but there are two reasons shrinking phase runs might stop: either the MT switched back to a growth phase, or the MT reached zero length.This latter possibility is incompatible with the idea that a Markovian switch rate should be the inverse of the average time spent in that state.Effectively, in our model, there is time spent in the "shrinking state" while stuck at length zero, causing a mismatch between the model and experimental definitions of the term "time spent in the shrinking state."As a result, we reserve the quantities λ ± s→g as degrees of freedom for the system, that we may use to satisfy other constraints.As we will see, special choices will be made to satisfy (C1).
By construction of the length-dependent catastrophe rate (Definition 2.1), constraint (C3) is immediately satisfied by selecting L 0 = L * and setting λ ± g→s = 1/τ ± g .Meanwhile, the selection L * = L addresses constraint (C2) if we make a complementary choice for the Michaelis-Menten constant F 1/2 .Indeed, the steady-state flux of tubulin into the microtubule pool can be expressed in two ways: first, G ± times the average growth rate of individual MTs that are in each growth phase, i.e. in terms of the experimentally-observable parameters v ± g ; and second, recalling Equation ( 6), in terms of the Michaelis-Menten uptake of tubulin, i.e. in terms of the steady-state value for available tubulin, F * = T tot − N L * and the maximum MT growth rates.
Together, we have the equation Now, since the ratios v + g /v +,max g and v − g /v −,max g are hypothesized to be the same (which is equivalent to assuming that the plus-and minus-end Michaelis-Menten constant is the same), we choose to separate the plus-and minus-end constraints and observe that v from which it follows (after imposing the target L * = L) that With this selection, (C2) will be satisfied.
While constraint (C1) may be the most obvious requirement to impose, it is the most difficult to address.By Proposition 3.3, we know that there is a unique steadystate solution, so if the number of degrees of freedom match the number of conditions implied at steady-state, then a solution should exist in principle.By virtue of the selection that L 0 = L * , the dimensional version of ( 13) implies that steady-state values for the proportion of MTs in growth phase at the plus-and minus-ends are Substituting these values in the left-hand side of ( 17) and setting them equal to the tubulin flux out of the MT pool, L * must satisfy If we separate the plus-and minus-end constraints, recall the definitions of ψ α in Equation ( 11) and L crit in Equation ( 10), and simplify, we have the system of equations This equation can be written more efficiently in terms of the dimensionless quantities Imposing the choice λ ± g→s = 1/τ ± g , we can rewrite this as Looking first at the plus-end version of this equation, δ + s can be expressed in terms of a generalization of the Lambert-W function: where W α (z) is the solution of the equation we −w α = z.
Using the right-hand side of (20) to find an expression for the exponential term, we can substitute into (19) to find a relationship between δ − s and δ + s .Namely, Since the velocity components of δ ± s are already prescribed, we can satisfy constraint (C1) by choosing L * = L and defining We have therefore determined all implied model parameters in Table 1 for our application to MT turnover in Drosophila dendrites.This completes our strategy for selecting model parameters for the stochastic model (see Figure 2), which we now numerically investigate in Section 4.

Comparison of the impact of rate-limited mechanisms
We investigate predictions of microtubule measurements based on stochastic simulation results in qualitatively different parameter regimes.As outlined in § 3.3, we wish to satisfy constraints (C1), (C2), and (C3) in order to have stochastic outputs match experimentally observable quantities (Figure 2).We set the total amount of tubulin to T tot = 1000µm and use the prescribed, experimentally-informed, and implied model parameters outlined in Table 1.The mean target and experimentally observable quantities (plus-end growth speed, plus-end growth duration, MT length) for T tot = 1000µm are shown with red dashed lines in the left panels of Figure 6 for the three values of the magnitudes of length-dependent catastrophe γ investigated in Figure 5: (a) γ = 0, (b) γ = 0.005, and (c) γ = 0.03.We show plus-end growth duration (1/λ * g→s,+ = 2 min) instead of the switching rate since it is a more interpretable quantity that can be obtained from experimental observations.The left panels of each parameter dashboard show that the distributions of the corresponding quantities that emerge from the stochastic model have averages that agree well with the target parameters (blue dotted lines).This provides evidence that the ODE model steady-state approximation provides useful insights for the stochastic model calibration.The right < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 L w w t 2 j r z N Y b 7 B R d dashboard panels show time-course information on tubulin allocation, minus-end position, and microtubule length.Note that each curve corresponds to the mean quantity across 10 simulations, and the cloud denotes the interquartile spread of each output.
The measurements used to generate the left panels of Figure 6 are extracted after 60 minutes of simulation time pooled over the 10 simulations, to ensure that the system has achieved steady-state.The top left panels pool all strictly positive speeds of the MT plus ends from the stochastic simulation, to be consistent with experiments that track the velocities when the plus ends are actively moving.Turning on the length-dependent catastrophe mechanism shows a tightening of the emergent speed distribution (see Figure 6b), due to the additional regulation of the microtubule lengths.The middle left panels show all the durations of growth at the plus ends of the microtubules, which are similar across the parameter ranges studied and show very good agreement with the target time spent in the growth state.The bottom left panels pool all strictly positive MT lengths.Interestingly, the length distribution for the γ = 0 case in Figure 6a predicts an exponential-like MT length distribution, while adding on the length-dependent catastrophe mechanism in Figures 6b and 6c predicts a length distribution with a peak around the desired target length.
The right panels of Figure 6 illustrate emergent properties of the stochastic simulations through time, as introduced in Section 2.1.The top right panels show the distribution of tubulin across available, unavailable, and microtubule pools.The length-dependent catastrophe mechanism limits the amount of tubulin in MTs and frees it up to the available pool, as expected.The middle right panels show the mean (red line) and the interquartile region (red region) of the MT minus end position trajectories, pooled across all MTs and all times.These panels show that, in the γ = 0 regime (corresponding to no length-dependent catastrophe), we predict a larger overall movement of the entire microtubules in the dendrite segment.Finally, the bottom right panel shows the mean (blue solid line) and the interquartile region (blue cloud) of microtubule length over time.Compared to the length distribution in the bottom left panels in Figures 6a-c, the microtubule lengths over time include zero-length microtubules.Therefore, the average length as a function of time will be slightly smaller than the average length in the corresponding histogram.As γ increases, the microtubule length distribution also tightens over time, where γ = 0 has the largest interquartile region.

The predicted response to tubulin manipulation
The previous section investigates how different length-regulating mechanisms impact the MT behavior in stochastic simulations with prescribed input parameters.The neuronal cell environment may however change through time and developmental stages, and in particular it may be characterized by different levels of tubulin.We are interested in predicting how these distinct tubulin levels may impact the MT growth and regulation.We assume that some growth characteristics of MTs would be fixed in varying tubulin environments.Namely, parameters F 1/2 , v± g , and τ ± g are likely innate characteristics of the cell that are unlikely to change as tubulin is varied.Therefore, we investigate a proposed tubulin titration experiment where we first parameterize our stochastic simulation at T tot = 1000µm based on the procedure in Section 3.3 < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 and observe how emergent MT properties are affected by subsequently varying tubulin.For each level of tubulin, we perform 10 simulations and aggregate the emergent distributions and time-course data as described in the previous section.
Figure 7 illustrates the impact of two extreme tubulin levels on the distributions of the emergent quantities and how the means of these distributions relate to the target parameters.Again, these distributions represent quantities taken after 60 minutes to ensure the system is at steady state.Microtubule length, plus-end growth speeds and plus-end growth durations distributions are shown for T tot = 700µm (light hue) and T tot = 4000µm (dark hue) with N = 20 microtubules.In the top panel, we show results for model simulations with no length-dependent catastrophe mechanism (γ = 0, blue).The emergent mean growth speed in the bottom left panel is below the target mean growth speed for the low tubulin level, and there is a wide spread of speeds observed in the simulation.We observe similar emergent distributions of the plus-end growth durations across tubulin levels, but also significantly different emergent length distributions for low and high tubulin.When tubulin is low, the emergent mean length is close to the target length, and most microtubules are shorter than the mean length.This is not the case for the high tubulin level, where the emergent mean length is much larger than the target length, and the exponential-like distribution shows much longer microtubules present in the system.
For small length-dependent catastrophe (γ = 0.005, red) in the middle row of Figure 7, the top right and center panels show that, for low tubulin, the target mean length and target plus-end growth duration are achieved, but the emergent growth velocity is below the target plus-end growth velocity.When tubulin is increased to T tot = 4000µm, the emergent growth velocity is larger than the target growth velocity and the distribution of speeds is much tighter.The emergent plus-end growth duration is also similar to the target duration, but we see slightly different emergent distributions of the plus-end growth duration.The emergent mean MT length in the top right panel is slightly larger than the target length and for both tubulin amounts, we find that the distribution of MT lengths is not exponential.When length-dependent catastrophe is increased further to γ = 0.03 (shown in green) in the bottom row of Figure 7, there is little difference between the low and high tubulin distributions, and the distributions of MT lengths are symmetric about the target mean.
We now aim to understand how the calibration and prediction observables change as we vary tubulin across a range of levels, so that T tot = {700, 800, . . ., 3900, 4000}µm in the stochastic simulation with N = 20 microtubules and L * = 35µm.We show the emergent mean and interquartile range of microtubule length in the left panel of Figure 8 for 10 simulations each, where the tubulin-limited only case is shown in blue, the γ = 0.005 case is shown in red, and the γ = 0.03 case is shown in green.For low tubulin levels (T tot ≤ 1000µm), the γ = 0 and γ = 0.005 simulations result in similar mean microtubule lengths.As tubulin increases, the emergent mean MT length continues to increase for γ = 0, while for γ = 0.005, the emergent mean MT length increases more slowly and appears to reach a plateau.The interquartile region for the γ = 0 case is much wider than the interquartile region for γ = 0.005, indicating that there is a larger variance of microtubule lengths for the case with no catastrophe.When γ = 0.03, length is very tightly regulated with a small interquartile region centered close to 35µm.Note that the mean MT lengths for γ = 0, γ = 0.005, and γ = 0.03 found in the stochastic simulations exhibit similar behavior to the analytical results shown in Figure 5 in Section 3.
Prediction observables also change as tubulin is varied.In the right panels of Figure 8, we show the average fluorescence turnover curve at each tubulin amount for (b) γ = 0, (c) γ = 0.005, and (d) γ = 0.03.The darker colors correspond to larger tubulin levels for each γ.As tubulin increases, the average fluoresence turnover curve moves up and to the right for both choices of γ.This indicates that more fluorescent tubulin is persisting in microtubules for more time at large tubulin amounts compared to smaller tubulin amounts.Comparing Figure 8b, Figure 8c, and Figure 8d, we see a larger spread in fluorescence turnover for γ = 0 compared to γ = 0.005 and γ = 0.03.As tubulin increases for γ = 0, it takes longer for fluorescent tubulin to turnover, indicating that for large amounts of available tubulin, the protein stays in the microtubules at the end of the simulation.For some tubulin levels, the timescale of turnover is consistent with previous experiments (Tao et al, 2016).This behavior is not seen for γ = 0.005, where all tubulin levels exhibit a decay in fluorescence that results in no fluorescent tubulin left in the activation window at 3 hours.When γ = 0.03, the turnover is even more tightly regulated for all levels of tubulin.A direct comparison with actual experimental data is a subject of further investigation.

Discussion
Uncovering minimal mechanisms of MT length regulation in living cells is crucial for ultimately answering questions about cytoskeleton dynamics and organization in healthy as well as post-injury settings.For example, in healthy neurons, a stably organized microtubule cytoskeleton needs to persist for the lifetime of the animal.Similar microtubule lengths should be present for years or even decades, while individual microtubules grow, shrink and turn over.We therefore develop a model that can be used in concert with available data on MT turnover dynamics from in vivo experiments.We consider a reduced differential equations model of MT dynamics at both ends using a traditional rate-balance compartmental model enhanced by terms that approximate intrinsically stochastic events like hitting the zero-length MT boundary.Steady-state analysis of this model informs our choices of parameters for a more extensive stochastic model of MT polymerization and depolymerization through switching between growth and shrinking states at both MT ends.The particular model parameterization we investigate here is motivated by experimental studies on MT turnover dynamics in segments of dendrites of Drosophila neurons (Feng et al, 2019;Rolls et al, 2021).
In contrast to prior studies, this modeling approach allows us to predict and match the spread of MT end speeds observed experimentally in the system inspiring this work.In addition, the model outputs predictions about the overall MT movement in the spatial domain, as well as about how tubulin protein may be allocated across MT and free tubulin pools in vivo.It is worth noting that our parameterization approach requires determining the rates of switching from shrinking to growth at MT ends based on the other available parameters.These rate estimates (and the associated run lengths spent shrinking) from both the deterministic and the stochastic models are different from those estimated experimentally.This could be due to the fact that run lengths may be cut off in experiments, in which data was captured for at most 20 minutes.This could either occur because of the MT switching back to a growth phase, or because the MT has reached zero length.
The most prominent contrast we draw in our modeling framework is between two physically plausible mechanisms of MT length control: tubulin-limitation and length-dependent catastrophe.We find that both mechanisms can yield MT length distributions with any given mean, but length-dependent catastrophe gives rise to much more variation in distribution shape.In the tubulin-limited regime, MT lengths appear essentially exponential in distribution, while in the length-dependent catastrophe regime, we observe anywhere from monotonically decreasing exponential-like distributions to unimodal Gaussian-looking distributions.Interestingly, this is consistent with experimental observations that have suggested that the steady-state length distribution of MTs is either exponential or Gaussian (Cassimeris et al, 1986;Schulze and Kirschner, 1986;Verde et al, 1990;Gliksman et al, 1992;Gardner et al, 2011Gardner et al, , 2013)).While the limited tubulin availability in neurons is not controversial, the exact molecular basis for the length-dependent catastrophe mechanism is less well established.Based on the modeling results presented here (in Sections 4.1,4.2),Gaussian-like MT length distributions with narrow peaks seem to suggest that both length-dependent catastrophe and tubulin limitations contribute to MT length control.Wider distributions would imply that tubulin availability may not be as large of a limitation to MT growth.The model provides a useful tool for helping distinguish between contributions of these mechanisms, and invites further study when experimental data on MT lengths in vivo becomes available.
Our stochastic model also predicts measurements that could be potentially observable in the future, such as the effective direction of movement of the MT ends.While the emergent dynamics of the minus ends requires further study, such insights on overall MT movement are typically not available from prior models.These insights can be further studied in the context of recent in vivo studies such as (Feng et al, 2019), which found that minus ends display sustained growth and that these dynamics are important for getting minus-end-out MTs into distal regions of healthy developed dendrites, as well as into developing and regenerating dendrites.
Similarly, the tubulin titration simulations provide a useful thought experiment for the impact of tubulin availability on MT regulation in the changing cell environment.In Section 4.2, we illustrate how low and high tubulin amounts influence the predicted growth speeds of MT ends as well as the MT lengths themselves.The two mechanisms of length regulation considered also generate distinguishable distributions of MT lengths and significant changes in the decay timescales for the model fluorescent turnover experiments.
The model thus generates testable predictions for experiments that manipulate tubulin in neurons, which have been proposed and could be used to validate hypothesized model mechanisms.In particular, if future experiments allow for control of tubulin levels in living cells, our model predicts that stabilization mean and variance of the MT lengths with increasing tubulin levels would suggest a significant contribution of the length-dependent catastrophe mechanism to MT length regulation.If photoconversion experiments are feasible at varying tubulin amounts, our model suggests that abundant tubulin regimes would be characterized by large changes in the timescale of fluorescent tubulin turnover with increasing tubulin levels (see results in Section 4.2).
A limitation of our current stochastic model is that it considers a fixed number of MTs, so that complete shrinking of a MT is followed by nucleation, or seeding, of a new MT.This modeling choice is motivated by the fact that nucleation is known to be heavily controlled in the cell.In particular, nucleation is heavily regulated in neurons (Chen et al, 2012;Hertzler et al, 2020), and plays a role in axon injury responses (Chen et al, 2012) and dendrite regeneration (Nye et al, 2020).However, it is not clear how robust this regulation is with respect to major changes in tubulin availability or other potential experimental manipulations.Depending on forthcoming experimental observations, allowing for a variable number of MTs could be an important future direction for enhancing our model.The proposed stochastic model and the parameterization procedure based on steady-state analysis of the reduced deterministic model can also extend to in vivo experiments in other animal systems, which are likely to be characterized by different MT lengths and polymerization speeds.units.To determine how small time step ∆t should be, we performed an experiment shown in Figure A1 to see the effect of time step size on average length.We see MT behavior is qualitatively the same for time step less than or equal to one second, so we utilize ∆t = 1 second.To determine how long each end spends in a certain state during time step ∆t, we first update the length-dependent catastrophe rate based on the length of the microtubules at the previous time step, t i , denoted as L i = L(t i ).Therefore, the catastrophe rate at t i+1 = t i + ∆t will take the following form where ϕ(x) = x − 1.Again, for small L and large γ, it is possible that λ min > λ ± g→s + γL 0 ( L(t) L0 − 1), so then the catastrophe rate will be λ min for short MTs.The updated catastrophe rates are then used to draw times spent in growth phase.
At each time step from t i to t i +∆t, each microtubule end is in either a shrinking or growth state.For each end of the microtubule, the time spent in growth and shrinking, denoted as t g and t s , respectively, are given by random variables drawn from the appropriate exponential distribution, where We then determine if the microtubule exits the time step in either growth or shrinking based on the idea that only one shrinking phase can occur in a time step ∆t (see Figure A2).For example, if a microtubule end entered in growth and t g > ∆t, then the microtubule would also exit the time step in growth (see Figure A2d).Similarily, if a microtubule end entered in growth and t g + t s < ∆t, then the microtubule also exits in growth, shown in Figure A2e.

Fig. 1 :
Fig.1: Overview of MT dynamics in a stochastic model simulation with two growth regulation mechanisms: tubulin-dependent growth and length-dependent catastrophe.Each microtubule has a length L and its plus and minus end are in either growth or shrinking states.In growth, the MTs utilize the available tubulin pool; while shrinking, they replenish tubulin into the unavailable pool.

Fig. 2 :
Fig. 2: Schematic outlining model parameters of the stochastic simulations, model observables that are used in parameterization of the stochastic model based on the deterministic model, and model outputs that provide predictions for experimental measurements.
Fig. 3: Outputs of a single stochastic simulation for (a) no length-regulating mechanism and for (b) length-dependent catastrophe (γ = 0.005) and tubulin limited growth (T tot = 1000µm) for N = 20 microtubules.The left panel shows positions of the plus end (blue) and minus end (yellow) of a sample MT throughout the entire simulation.Black dots indicate times when the MT reaches length zero.The right panel shows the mean and interquartile range of the MT lengths (orange) and each MT length (grey) throughout the simulation time.The horizonal line denotes the target MT length, L * = 35µm.Note the difference in the y-axis range for both panels in (a) and (b).

Fig. 4 :
Fig. 4: Cartoon of photonvertible tubulin experiment, inspired from (Rolls et al, 2021).(Top) Tubulin in existing MTs in a window of interest is fluoresced at the photoconversion time.(Bottom) The intensity of the remaining fluoresced tubulin in the same window is tracked at subsequent times.

Fig. 5 :
Fig. 5: (Left) Various evaluations of the functions H 0 (ℓ) and H 1 (ℓ), see Equation (16).We use the baseline parameter set in both cases, with the modifications that H 0 is evaluated over a range of γ (blue, red, and green), while the function H 1 is evaluated over a range of total tubulin (light gray to black as T tot increases.)The steady state values (filled circles) occur at intersections of H 0 and H 1 .(Right) The steady state L * is displayed as a function of T tot for the fixed values of γ.The displayed steady state values (filled circles) correspond to those displayed on the left.

Fig. 7 :
Fig. 7: Distributions of target quantities for tubulin titration experiment with tubulin amounts T tot = 700µm (light hues) and T tot = 4000µm (dark hues) for N = 20 MTs with input parameters set for γ = 0 (top, blue), γ = 0.005 (middle, red) and γ = 0.03 (bottom, green).Target and emergent mean values for various tubulin levels are denoted as in the dashboard figures.

Fig. 8 :
Fig. 8: A contrast of responses to manipulating the level of tubulin in the cell for N = 20 MTs when length-dependent switching to shrinking is present (red curves, γ = 0.005, green γ = 0.03) and is not (blue curves, γ = 0).(Left) Mean and interquartile region of average microtubule length as a function of total tubulin (µm).(Right) Average relative fluorescence intensity in models of tubulin photoconversion for various total tubulin levels.Colors in (b), (c), and (d) correspond to γ = 0, γ = 0.005, and γ = 0.03, respectively.
Fig. A1: Average MT length over time for different stochastic time steps, ∆t for γ = 0, L * = 35µm, and N = 20 MTs for one realization.

Table 1 :
Parameters that appear in the model as well as experimentally observed quantities that inform our choices.Parameters marked with ‡ are preliminary estimates from the Rolls Lab.A full description of our parameterization strategy is addressed in Section 3.3.